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ABSTRACT 

High redshift sources suffer from magnification or demagnification due to weafc gravi- 
tational lensing by large scale structure. One consequence of this is that the distance- 
redshift relation, in wide use for cosmological tests, suffers lensing-induced scatter 
which can be quantified by the magnification probability distribution. Predicting this 
distribution generally requires a method for ray-tracing through cosmological N-body 
simulations. However, standard methods tend to apply the multiple thin-lens approx- 
imation. In an effort to quantify the accuracy of these methods, we develop an in- 
novative code that performs ray-tracing without the use of this approximation. The 
efficiency and accuracy of this computationally challenging approach can be improved 
by careful choices of numerical parameters; therefore, the results are analysed for the 
behaviour of the ray-tracing code in the vicinity of Schwarzschild and Navarro- Frenk- 
White lenses. Preliminary comparisons are drawn with the multiple lens-plane ray- 
bundle method in the context of cosmological mass distributions for a source redshift 
of Zg = 0.5. 



Key words: cosmology: theory 
scale structure of Universe 



gravitational lensing - methods: numerical - large 



1 INTRODUCTION 

Weak gravitational lensing of distant galaxies by large-scale 
structure leads to the distortion of their images. In an inho- 
mogeneous universe, each line of sight probes a slightly dif- 
ferent integrated mass and is sheared differently; overdense 
regions result in magnification and underdense regions cause 
demagnification. The result is that for a given redshift z, a 
source is magnified (or demagnified) by some amount /i, rel- 
ative to the case of a perfectly homogeneous universe on all 
scales. This magnification has an associated probability dis- 
tribution that reflects the existence of structure on a range 
of scales. 

Large surveys for galaxies, quasars and supernovae suf- 
fer a two- fold magnification bias as a result of the same phe- 
nomenon. When a source is magnified, its total surface area 
on the sky is increased; consequently, for a given area of 
sky observed, the region of the source plane being sampled 
is decreased. On the other hand, magnification will push 
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otherwise too-faint sources above the observa tional thresh- 
old i n flux-limited surveys (Canizares 198ll : [Turner et al.l 
Il984 ): this is p articularly imp o rtant for optically-selected 
quasar surveys ([Peacock Il982l : iBartelmann fc Loebl [l99^ : 
IWvithe et al. 12011 ). Together, these two effects have coun- 
teracting influences on source counts, but generally do not 
cancel each other out. Additionally, for sources that are 
demagnified the effect is reversed; indeed, the majority of 
sources are de-magnified relative to the case of a perfectly 
homogeneous universe. The net effect on source counts, ob- 
served luminosity functions, source redshift distributions 
and any resulting bias depends on selection procedures, in- 
trinsic source properties and the probabilistic lensing effec t 
(|Dver fc Roedeilll974l : |Peacocklll982l : iLe Brun et all I2OO0I ') . 
The intrinsic source luminosity function and magnification 
probability distribution function (hereafter ^PDF) are un- 
knowns; the former is ge nerally modelled as a Schechter 
function l|Schechteij Il976l ). An appropriate model for the 
latter — which may be produced by ray-tracing through 
N-body simu lations — is the subj ect of the present study 
(see Sec.l ^ . lHamana et all (|20od ) and lJain fc Limal (|201ll ') 
have demonstrated that a power-law tail in the /iPDF signif- 
icantly changes the shape of the bright end of the luminosity 
function and generates a considerable magnification bias for 
high-redshift sources. 
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Gravitational lensing causes magnification or demagni- 
fication of Type la supernovae (SNla), resulting in a scat- 
ter in the inferred distance-redshif t distribution, p articu- 
larly for high-redshift sources (jJonsson et al.l 12008 ). The 
lensing-induced scatter introduces a bias and uncertainty 
in the inferred values of the matter density parameter, 
f^M.o, the deceleratio n parameter, go and the dark energy 
equa t ion of state, w ( Wa,mbsganss et al] Il997l : lAstier et al.l 
I2OO6I : IJonsson et al.l [20081 ). For example, in the case of the 
z = 1.7 supernova SN1997ff, the estimated magnification 
would easily bridge the gap between evidence supporting 
a ACDM universe and that supporting an empty universe. 
The lensing effect can be 'averaged out', for supernovae sam- 
ples, by combining measurements from a sufficiently large 
number of sources a t similar redshifts ||Ho1z fc Wald|[l998l : 
iHolz fc Lindeill2005l ). However, a similar scatter will be in- 
duced by lensing in the Hubble diagram constructed from 
measu rements of gravi t ational waves, the so- c alled standard 
sirens l|Markoviclll993 : IHoIz fc Hugiie3l2005l : [ Jonsson et al.l 
l2007h . Efforts are now underway to provide an appropriate 
method for correcting the observed brightnesses of individ- 
ual objects for magnification (e.g. [Gunnarsson et al. 2006; 
IJonsson et al.ll2009l : ishapiro et al.ll2010l : lHilbert et ahlboill ^. 



By accounting for the non-Gaussian shape of the /iPDF, one 
is best able to correct for the effects of len sing statistically 
l|Hirata et al.ll2010l : [Shang fc Haimanll201ll ). 

There are many approaches to determining the /iPDF 
for a given cosm ological model. The optical scalar equa- 
tions l|SachJl96ll ). which describe the evolution of the cross- 
section of a small beam of light, have been applied to var- 
ious mass distributions leading to useful redshift-distance 
relations in the limiting c ase when the line of sight is far 
away for inhomog eneities l|Kantowskilll969l : iDver fc Roederl 
1 1974 iDveil [ 19771 ). The equations have also been applied 
to an infinitesimal beam transported through a general- 
ized inhomogeneous universe that, on average, satisfies the 
Friedmann-Lemaitre-Robertson- Walker (FLRW) geometry 
in ord er to re-derive th e maximal angular diameter dis- 
tance l|Seitz et al.l Il994l ). A simple integration of the op- 
tical scalar equations would be possible in the case of a 
known metric, but that is the crux of the problem - the 
local metric for an inhomogeneous universe has no general 
(analytic) solution. Crude model universes incorporate inho- 
mogeneities^ and descr ibe their effects o n light pro pagation 
(e.g. iDver fc Roeder.igT^ : lLinder|[l988l : iFutamasc fc Sasakil 
Il989l ). However, when the cosmological structure being 
probed is highly non-linear and lines of sight have small im- 
pact parameters, studying light deflection requires the use 
of numerical techniques. 

The non-linear and hierarchical growth of the large- 
scale structure is generally modelled by cosmological N-body 
simulations, with the propagation of light and its subse- 
quent lensing computed by ray-tracing through these sim- 
ulations. Some methods consider only one dominant lens 
in each line of sight; the lensing object is thin compared 
to the distances between observer, lens and source - this 
is known as the thin-lens approximation. However, multiple 
lenses coincidental along the line of sight may be responsible 
for a lensing event jWambsganss et al]|2005l : [Das fc Ostrikej 
I2OO6I ). The existe nce of multiple lens p lane also accounts 
for the findings of iPremadi et al.l l|200ll ): that the magnifi- 
cation probability, P(/i > 1) is mostly independent of the 



parameter that qua ntifies the normalis ation of the matter 
power spectrum, erg . iRauch et all l|l992l ) identified a feature 
of the /iPDF: a 'bump' that was only evident when con- 
sidering two-dimensional projections of matter; they specu- 
lated that changes to the caustic structure resulting from the 
projection was somewhat responsible. Th e presence of thi s 
caustic- induced bump was confirmed bv iLee et al.l |l993), 
who derived a semi-analytic expression for this feature in 
the limit of low-optical depth. Therefore, in the context of 
cosmological structure, the thin-lens approximation should 
be replaced, at the very least, by the multiple-lens plane 
approach (see Sec. 13. ip . where large volumes of matter are 
proje c ted onto a seri es of lens planes (|Blandford fc NaravanI 
1 19861 : iKoyneil 1 19871) . The ra y-shooting method , deve l- 
oped bv iKavser et al.l ll 19861): ISchneider fc Weisi l|l988t ): 
IWambsgansd (|l990l ): IWambsganss et al.l (| 19981 ). embodies 
this approach. Subsequent introduction of tree-methods to 
measure the defiection angle at each lens-plane have pro- 
duced efficient algorithms. Early ray-shooting techniques 
were applied to a random distribution of point masses, 
which therefore did not incorporate the intricacies of 
clustering properties wit hin various cosmological models 
JSchneider fc Weisd 1 19881 : iPaczvnski fc Wambsgaiisi Il989l : 
iLee fc Paczvnskilll990l )~ Later studies coupled the multiple 
lens-plane approach with mass distributions taken from N- 
body simulations to study the effect and magnitude of grav- 
itational lensing (e.g. Jaro szvnski et al. 1990). 

While it is tempting to continue using the previous 
methods because of their simplicity, one must first quantify 
the effect that the various approximations have on the result- 
ing predictions. Galaxy redshift surveys su ch as the CfA sur - 
vey l|Davis et al.lll98a ) and the 2dFGRS (|Cole et al.ll2005l ) 
have revealed that elongated structures exist in the form of 
filaments, which stretch across large voids between galaxy 
clusters. These are also evident in cosmological simulations. 
If such filaments are projected onto lens planes, the result- 
ing magnification would be overestimated. Though earlier 
computational limitations, such as memory and processing, 
warranted the need for these simplifications, we are now en- 
tering an era where more accurate methods are within our 
reach. 

Only a few studies have numerically integrated the 
null geodesic equations from observer to source. The ear- 
liest of these studies assumed metrics that were approxi- 
mated via a simplified model for inhomogeneities and de- 
rived distance-redshift relations that were compared to 
the Dyer-Roeder approximations (Futamase fc S asakil 1 19891 : 
Watanabc fc Tomita 1990). Kasai et al. (1990), drew at- 
tention to the spread in angular diameter distances given 
a small enough beam-size and the fact that the average 
value was consist ently lower th an the solution for a homo- 
geneous universe. iTomital l| 19981 ) numerically integrated the 
null geodesic equations through N-body simulations, albeit 
at low-resolution, finding that various cos mologies exhibited 
differe nces in angular diameter distances. ICouchman et all 
(| 19991 ) developed a three-dimensional algorithm for mod- 
elling weak gravitational lensing; comparing the two- and 
three-dimensional shear, they found that the projection of 
structure on 100 h~^Mpc scales led to errors of up to 9 per 
cent d epending on the redshift of the lens box. I Vale fc Whitg 
(|2003| ) performed three-dimensional ray-tracing by comput- 
ing the defiection angle for rays many times along their path. 
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advancing them in step-sizes of Lso^/Ng along the line of 
sight. In their study, they used a relatively small number 
of grid points along the line-of-sight Ng = 32, such that 
the step-size was approximately 10 h~^Mpc, but boxes of 
half that size were used for numerical tests. They found 
that their convergen ce maps were unchanged to 0.1 per cent. 
IWhite fc Hul (|2000l ) run simulations of structure formation 
while simultaneously propagating photons through said sim- 
ulation; they were able to produce convergence fields from 
the light rays they follow. 

What has been missing thus far is a one-to-one com- 
parison of results that take a multiple-lens approach and a 
direct numerical integration of the null geodesic equations; 
this comparison is carried out in the present work. In Sec- 
tion [21 the basic theory and notation of gravitational lensing 
and the relevant statistics are presented. We introduce the 
method of numerically integrating the null geodesic equa- 
tions in Section |3] with specific reference to its multiple lens 
plane counterpart; in Section IH we establish the accuracy 
and limitations of the method. We describe the simulations 
used for modelling the cosmological mass distribution and 
compare statistical predictions determined by the two meth- 
ods in Section [5] Finally, a summary of our findings is pre- 
sented in Section [6] 



2 GRAVITATIONAL LENSING 

Here, we introduce the conventional notation for the com- 
ponents that describe the effect of gravitational lensing on 
the shape of a beam to first order. Each background source 
experiences lensing in the form of convergence, k, and shear, 
7. Convergence is the isotropic Ricci focusing of a beam due 
to enclosed matter, while shear is the tidal stretching of the 
beam along a particular axis due to asymmetric matter dis- 
tribution; together, they serve to increase the area of the 
source on the sky, resulting in its magnification, due to the 
conservation of surface brightness. The two-dimensional ef- 
fective lensing potential, '4>2d, is given by 



'^2d{'x) ~ — / d'^2;'K(x')ln |x — x'l 

71" Jk2 



(1) 



Here x is a dimensionless vector formed by scaling the image 
position, ^: 



Co' 



(2) 



For a single thin lens, k is equivalent to the scaled surface 
density of the lens, and is related to the gravitational po- 
tential via the Poisson Equation 



1^2 , S 

-V V2D = K = , 



(3) 



where the critical surface density for gravitational lensing is 
given by 



Ds 



(4) 



Thus lensing is characterised by Ds , Dd and Dds , the angular 
diameter distances from the observer to the source, from the 
observer to the (thin) lens, and from the lens to the source 
respectively. Total shear can be written in complex notation; 



its dependency on the two orthogonal components are given 
by: 



7 = 71 + ^72, 



(5) 



where the two components are linearly related to the second 
derivatives of the (projected) gravitational potential along 
two orthogonal directions by 



71 



- ^,22) 



and 



72 = tp.-i 



(6) 



where the indices after the comma denote partial differenta- 
tion and we temporarily drop the subscript 2D for clarity. 
The deformation of the beam is described as a mapping from 
the source plane to the image (observed) plane. The Jaco- 
bian. A, of the lens mapping is a real and symmetric 2x2 
matrix given by 



1 - K - 71 
— 72 



-72 
1 — K + 71 



(7) 



The fiux magnification of an image is given by the inverse 
of the determinant of the Jacobian, so 



(1 



l7l = 



(8) 



One may now consider the effect of lensing on the appar- 
ent position of the image of the source of light. The scaled 
deflection angle is the gradient of the lensing potential: 



020 = VV'215 
1 



— — d ^'^(x') 



(9) 
(10) 



The gravitational lens equation is therefore given by 

l3 = 0-a2D{0), (11) 

where (3 is the angular source position and is the angular 
position of the image on the sky. The deflection angle, a, is 
related to its scaled counterpart by: 



a2D = 



D.Dd 



■OL-2U- 



(12) 



2.1 The magnification probability distribution 

The probability that a source at redshift z would be magni- 
fied by an amount within the interval [/i, /i-f d/i] is p(/i, z)d[i,. 
It satisfies : 



(13) 



When the /iPDF is convolved with intrinsic luminosity dis- 
tributions for standard candles, the result describ es the ob- 
serve d spread in magnitudes. Flux conservation l|Weinberel 
1 19761 ) demands that p(/i, z) satisfies 



tJ.p{fi, z)d^ = 1. 



(14) 



There exists a minimum magnification, indeed a minimum 
convergence, which corresponds to a line of sight that en- 
counters no matter between obs erver and source; this is also 
referred to as an empty beam (|Dver fc Roedeilll972h . The 
distribution function peaks at values below ji = 1 and is 
highly skewed towards high magnifications jHamana et al.l 
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I2000I 1. For low redshift sources and/or small lensing opti- 
cal depths, the high magnification tail will exhibit a power 
law trend p(^) oc which results in a formally divergent 
standard deviation ( Vietri fc Ostr ikcr 1983). This result was 
originally derived analytically for a random distribution of 
compact l enses where only one lens dominates; however, 
[Pdl (|l993l ) have noted that for higher optical depths, the 
slope of tail may become shallower. The precise shape of the 
distrib ution depends o n the assumed density profile of the 
lenses ||Yoo et al.ll2008l '). Th e spread in this dis tribution in- 
creases with source redshift (|Babul fc Leelll99ll ). For a fixed 
source redshift, the shape of the distribution, particularly 
for low magnifications, depends mostly on the cosmologi- 
cal parameters as an d the vacuum density parameter, JIa.o 
l|Premadi et al.ll200ll '). The high magnification region suffers 
from low number statistics and since it represents the effects 
of strong lensing, cannot be probed by ray-bundle methods, 
which we now discuss. 



3 RAY TRACING 

3.1 Multiple lens plane methods 

The multiple lens-plane method requires the total three- 
dimensional mass distribution to be sliced up into contigu- 
ous boxes; each box is then projected onto a plane perpen- 
dicular to the line-of-sight, usually placed at the centre of 
the box. In backward ray-tracing methods, light rays are 
propagated from the observer to the source with deflections 
only occurring in successive lens planes. The deflection at a 
given plane is the result of matter within that plane only. 
The formalism for the resulting deflection can be written in 
terms of t he m ultiple Icns-planc equatio n, derived and de- 
velope d bv lBlandford fc NaravanI (Il986l ') , ISchneider fc Weiss! 
l|l988ll and lJaroszvnski et al.1 (|l99(]| ): the source position, r], 
after deflection by A'^ lens planes: 



Do 



(15) 



where is the position of the ray in the i"^ lens- plane; ^1, 
therefore, is the position of the image. Dos, Doi and Dis 
are the angular diameter distances from the observer to the 
source, from the observer to the flrst lens-plane, and from the 
i*^ lens to the source respectively. Notice that high density 
regions in a plane may not necessarily be dense in real space. 

One approach to studying the statistics of gravitational 
lensing is the ra y-bundle me thod fRBM) , developed by 
iFluke et al.l (|l999l i. and see also lFlukd (| 19991 ). Here, the mul- 
tiple lens-plane method is used to propagate light rays from 
the observer to a given source plane. Instead of employing a 
grid-based technique for calculating magnifications across 
the source plane, the app roach favou red b y other back- 
ward ray-tracing codes (e.g. jjain et al.ll2000l : [Premadi et al.l 
the RBM models each individual line of sight as an 
'infinitesimal bundle' of light rays around a central ray and 
follows this bundle back to the source plane. An initially 
circular image is distorted by the time it reaches the source 
plane as a result of convergence and shear. These quanti- 
ties can be determined for this specific line of sight as the 
image-source association is maintained. In the RBM, N-rays 
in a regular polygon represent a circular image and their 



positions at the source plane are fitted with an ellipse; the 
Jacobian matrix (Eqn. [7)l is determined for each bundle and 
solved to determine fi, k and 7. 

The ray-tracing method presented in this work uses the 
RBM design of a ray-bundle with 8 rays, but does away 
with the multiple lens-plane treatment of the lensing mass. 
Instead, the evolution of the cross section of the bundle is de- 
termined by integrating the null geodesic using a numerical 
gravitational potential obtained from N-body simulations. 

3.2 The weak field metric 

The Friedmann-Lemaitre- Robertson- Walker (FLRW) geom- 
etry describes a universe that is homogeneous on all scales 
The FLRW line-element for a fiat geometry is: 



ds 



+ w [t) [dx + X (de'^ + sine < 



(16) 



where S and (p are comoving coordinates. The metric al- 
lows an evolving scale factor, R{t). From here onwards, we 
quantify the scale factor relative to its current value, i.e. at 
t = 0: 



a{t): 



m 

RiO)- 



(17) 



We assume that the inhomogeneities present in large scale 
structure are small enough to be represented as a perturba- 
tion to the background FLRW metric, and do not falsify the 
large-scale predictions made under FLRW geometry. The 
resulting line element is: 

ds^ = -[1 + 2ip{t,x)]dt^ + a^(t)[l - 2ilj{t,x)]dx'^ (18) 

where i denotes the three spatial dimensions. The weak-field 
metric is applicable where tp{t, x) <^ a{t). The gravitational 
potential, which is defined with respect to the local per- 
turbation from a smooth background (see Sec. 13. 3p . can be 
decoupled 



^lj(t,x) 



{x) 



a{t) 



(19) 



such that tp is defined in physical units and in co-moving. 
The Christoffel symbols for this metric, presented in Ap- 
pendix [A] are dependent, not only on the gradients of the 
gravitational potential <p, but itself, which is defined based 
on the mass perturbation from a smooth background. The 
Geodesic Equations (see Eqn. [27] below), which are the sec- 
ond order differential equations for the four coordinates, are 
then constructed. 



3.3 The gravitational potential 

The perturbation field: 

p{t,x) - p{t) 



S{t,x) 



(20) 



relates the local density, p{t, x), to the mean matter density, 
p{t), the latter of which is given by 



p{t) = QMPc{t). 

The critical density, pc(t), is given by 

3^2 (t) 



Pc{t) = 



87rG 



(21) 



(22) 
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where H{t) is the Hubble parameter: 

m 



H{t) 



(23) 



The perturbation is related to the gravitational potential 
ip(t,x), by 

^it,.)^-'-^pit)f d.-^. (24) 



The derivatives of the gravitational potential are then 

dth(t,x) AttG , , f , n^, , x — x' 
ax^ ,b3 X 



(25) 



One can relate the numerically determined (j) to the grav- 
itational potential via Eqn. [19] and derivatives in a similar 
fashion: 



dtp{t, x) 
dx^ 



1 d(l>{x) 
a{t) dx^ 



(26) 



Constructing the grids that represent the perturbation field 
and its derivative requires a few steps. Firstly, if the lensing 
mass is discretized into particles, then the mass within a 
cube of co-moving side length Lbox is assigned to the nodes 
of a regular g rid using the Cloud-in-Cell algorithm (CIC; 
iHocknev fc Eastwood 1988) in three dimensions; the mean 
matter density is then subtracted off. A Fast Fourier Trans- 
form is applied to convolve the mass distribution with the 
appropriate kernels to determine the gravitational field and 
its derivatives according to Eqns. [24l and l25l using the popu- 
lar software package 'Fastest Fourier Transform in the West' 
(FFTV0). Depending on the lens, a periodic mass distribu- 
tion may be implied, or not. If not, the density grid is zero- 
padded before performing the FFT convolution. At many 
stages during the integration, the local values of the field 
and its derivatives needs to be evaluated, so an interpola- 
tion scheme was written for this purpose; it is described in 
detail in Appendix iBl 

There are two numerically intensive parts of the ray- 
tracing method: the first is the set of Fast Fourier Trans- 
forms required to calculate the gravitational potential and 
its derivatives at every grid point on a 3D mesh; the second 
is the interpolation required at each time-step to determine 
the values of the same quantities at the exact position of the 
light ray. 

Note also, that the Fourier-grid resolution sets a lower 
limit to the scale of structure probed, as the mass is 
smoothed over this grid scale. This is assuming that this 
scale is reasonably larger than the softening length employed 
in the N-body simulations used to model the cosmological 
mass distribution. 



3.4 Three-dimensional ray-tracing 

Our approach is based on avoiding the approximation of the 
multiple lens plane method. Instead, the path of a photon 
is numerically integrated from the Geodesic Equation: 



dA2 



^ http://www.fFtw.org/ 



= -F' 



dx^d3£_ 

IxHx 



(27) 



Here F^^ denotes the Christoffel symbols associated with the 
metric, x'' represents any of the four coordinates specified by 
the superscript i, and A is the affine parameter. 

The four second-order differential equations are reduced 
to eight coupled first-order differential equations, which are 
integrated with the affine parameter A as the dependent 
variable. A classical fourth order Runge-Kutta integration 
scheme with fixed step-size was written and used to perform 
the integration. The subtlety is that it is the -coordinate, 
rather than the affine parameter, that determines the bound- 
aries of the integration. The exact evolution of the affine 
parameter will be different for each ray-bundle, and is not 
known in advance. Therefore, the integral is repeated in 
small steps of A, AArk, until the x^ reaches the value re- 
quired at the source plane Dc{zs)- The step-size is chosen 
such that the estimated resultant step in x'^ is a fixed frac- 
tion, /rk, of the Fourier-grid resolution: 



AAe 



fnKLbox 



ufptx: 



(28) 







where the overdot denotes the derivative with respect to A, 
the extra subscript denotes values at t = 0, and wpft 
denotes the number of points along one side of the Fourier- 
grid (see Sec. 13. 3p . If the ray overshoots the source plane, a 
simple linear interpolation, using the current and previous 
positions, is used to find the final source position. The use of 
a fixed step-size is deemed appropriate as the values of the 
Christoffel symbols are interpolated from gridded values, as 
described in Sec. 13.31 There would be little information gain 
from step-sizes much smaller than the grid resolution. Never- 
theless, the effect of the choice of step-size is analysed, along 
with other numerical parameters, with results presented in 
Sections |4T1 and IT2] 



3.5 Evolution of the scale factor 

For certain cosmologies, the Friedmann Equations can be 
solved to find the specific time dependence of the scale 
factor; for a spatially fiat, radiation-free Lemaitre model 
(fik ~ 0, flu + ~ 1), with JIa > 0, this dependence 
is given by : 



a{t) 



^M,0 . , 2 

smh 



1/3 



where to is given by 



to = 



sinh 



'm,o 



(29) 



(30) 



Eqns. 



were derived from Eqn. 15.36 in 
the solution for a spatially fiat, matter- 



iHobson etHI |20& 
only Lemaitre model, but with t used to denote lookback time 
instead. 

As the photon traverses a single lens box, the scale fac- 
tor will evolve (i.e. decrease) although the co-moving scale 
of the structure does not change appreciably over this time- 
scale. The value of the scale factor is required to evaluate 
the Christoffel symbols, so we evolve the lookback time for 
the photon as well as the spatial coordinates and use this to 
determine the scale factor at each position throughout the 
box. 
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4 COMPACT LENS MODELS 

Here, we present results of ray-tracing in the vicinity of sim- 
ple lenses. Both of the lens models used are fairly compact, 
and so the thin lens approximation applied in the analytic 
solution is suitable. The analytic solutions should be read 
as magnification for the image at that location, rather than 
the total magnification of the source. This is an important 
distinction. The ray-bundle approach only follows the con- 
gruence corresponding to a single image, and therefore can- 
not account for the total magnification of multiply imaged 
sources; this is one of the reasons why it should not be used 
to model strong lensing. Yet for these simple lenses, we are 
still able to test our results against the known magnification 
of a single image for an otherwise strongly lensed source. 



4.1 Schwarzschild lenses 

The Schwarzschild lens, a singular density point, is the sim- 
plest lens one may study. The analytic solution for the mag- 
nification for an image at any location is well known and is 
presented in Appendix [C] Since the gravitational potential 
here is equivalent to a kernel (see Eqn. I24|l multiplied by the 
mass of the lens, we directly compute the gravitational field 
at each grid point, before performing the ray-tracing. The 
solution for this lens is scale-invariant, so by increasing the 
mass of the lens, we may essentially increase the resolution 
of the grid, and make it possible to study the behaviour of 
ray-bundles that pass near (or inside) the Einstein radius. 
In cosmological simulations, the RBM explicitly avoids this 
as the images near the Einstein ring are part of multiple set 
of images, and the images inside the ring are the demag- 
nified images which contribute only a small portion of the 
magnification of the associated source. However, for testing 
purposes, we allow lines of sight that approach the Einstein 
radius since the analytic solution is known for each individ- 
ual image. In the fiducial test case, the Schwarzschild lens 
has a mass of lO^^M©. It is placed at a redshift of zl = 0.35, 
giving it an Einstein radius of 65'.' 6 for a source redshift of 
Za = 0.8, which corresponds to ~ 12 mesh-points on the 
256^ Fourier-grid. 

The ray-tracing code developed here is able to repro- 
duce the desired lensing distortion very well. The plots on 
the left hand side of Fig. [T] show that even high magnifica- 
tion events, which occur when the image is near the Einstein 
radius, are recovered by the numerical method. The numeri- 
cal error, shown on the right hand side, increases for images 
that are closer to this radius where the higher-order lensing 
effects, like flexion, are expected to play a role and an eight- 
ray bundle has limited application regardless of the choices 
of other numerical parameters. Note that in each of the tests 
shown in the Figures[l]-|4l the results break down in the very 
centre of the lens. We choose 3 per cent as an error thresh- 
old to mark the sharp upturn in the average error at small 
impact parameters. The error surpasses 3 per cent when the 
impact parameter is smaller than a limiting radius, hereafter 
referred to as the 'minimum reliable radius', or MRR. The 
MRR can be used as a measure of numerical accuracy and 
the effect of numerical parameters. For example, the top row 
of Fig. [T] shows the result of decreasing the number of inter- 
polation points so that just one Fourier-grid node on either 
side of the current position is used to find the local values 



of the gravitational field and its derivatives. In this case, 
the MRR expands out to a larger impact parameter where 
only low magnification regions (^ < 1.5) are recovered with 
< 3 per cent accuracy. However, once the number of inter- 
polation points is increased beyond this, this parameter has 
little influence on the accuracy of the results. Likewise, on 
the bottom row of Fig. [T] we show the result of reducing the 
Runge-Kutta step-size; a step that is approximately half the 
size of the Fourier-grid resolution (/rk = 0.5) is sufficient 
for accuracy at moderate magniflcations (1.5 M ^ 6). De- 
creasing this parameter has negligible effect on the results. 

4.2 NFW lenses 



Vario us studies of cosmological simulations (|Navarro et al.l 
I1995') have found that dark matter haloes on galactic and 
cluster scales have mass distributions that are well described 
by a Navarro-Frenk- White (NFW) proflle: 



p(r) = 



5cpc 



(r/r,)(l + r/r,)2' 



(31) 



Here, pc is the critical density (see Eqn. I22|l at the halo 
redshift; r^oo is the radius within which the mean density 
of the halo is 200pc; '"s = '"200 /ca is the characteristic scale 
radius, which marks the transition in the the slope of the 
profile; and ca is the dimensionless concentration parameter. 
Finally, the characteristic overdensity is given by: 



5c 



200 



4 



3 ln(l+CA)-CA/(l + CA)' 



(32) 



We use this profile to compare how the multiple-lens plane 
approach and the three-dimensional approach model lensing 
around galaxy and cluster haloes, where most of the high 
magnification events will occur. 

Lenses with the NFW profile were modelled by discretis- 
ing the mass contained in one virial radius into a number of 
particles of equal masses and constructing a fake simulation 
output in the GADGET format. The numerically determined 
gravitational potential is therefore not equivalent to the an- 
alytic solution, but suffers from discretisation effects just as 
a simulated halo would. The spherically symmetric lens is 
divided into a fixed number of radial bins of equal width, 
extending from the centre to the virial radius. Each particle 
is randomly placed in one of the radial bins with a probabil- 
ity proportional to the mass within that bin; the bin mass is 
found by appropriately integrating the density profile, given 
in Eqn. 1311 The lenses modelled for this purpose have a virial 
mass, and total mass, of Mvir/MQ = 10^'' and concentra- 
tion parameter ca = 7.2. At a (lens) redshift of zl = 0.35 
and a virial overdensity of Ac — 200, their virial radii are 
0.55 h~^Mpc. The results presented here are for a source 
redshift of zs = 0.8. Figures [2] - 14] show the magnifications 
determined by ray-tracing over a range of impact parame- 
ters changing one numerical parameter relative to a fiducial 
choice. The ray-bundles are allowed to encounter the lens at 
a distance of up to twice the virial radius, which corresponds 
to an angular separation of 343" . The analytic solutions for 
image magnification and shear by a lens with the NFW pro- 
file are presented in Appendix [D] Since the model lens is 
truncated at one virial radius, the metric on the exterior is 
given by the Schwarzschild metric, by Birkhoff's theorem. 
Thus the analytic solution shown for p{r) switches to that 
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Figure 1. The radial dependence of magnification {/i; left) and the percentage error in the numerical magnification {Afi/ii; right) when 
ray-tracing past a Schwarzschild lens. Different values of numerical parameters are compared; results are shown for (Top) Number of 
points used to interpolate the local potential field and its derivatives: ni„i = 1 (blue squares), ni„t = 2 (red circles) and riint = 3 (pink 
triangles); (Bottom) Runge-Kutta step-size that is: 1/2 the Fourier-grid resolution (red circles) or 1/3 the grid resolution (blue squares). 
In each case, 10^ lines of sight have been distributed and binned uniformly with log(r/rvir) and la error-bars are shown. The analytic 
solution (black dots) at each impact parameter sampled is included in the panels on the left. In the panels on the right, the dotted 
horizontal line marks the 3 per cent error threshold. 



appropriate for a Schwarzschild lens outside this impact pa- 
rameter. 

Between the MRR and the truncation radius, the error 
remains fairly constant at 1 per cent. Outside the truncation 
radius, the error drops sharply by an order of magnitude, 
which imphes that the removal of mass outside this radius 
is primarily responsible for the aforementioned 1 per cent 
error. The size of the ray-bundle, or equivalently, the radius 
of the image or beam, has negligible effect on the ray-tracing 
results (see top row of Fig. [2]). By varying the grid size (see 
top row of Fig. [3]), we notice that the MRR has a very strong 
dependence on the Fourier-grid resolution. Our interpreta- 
tion is that the interpolation scheme used to measure the 
local gravitational field would smooth over the central cusp, 
causing the large errors here. We conclude that if small-scale 
structure is responsible for a caustic (high magnification) 
then this will be washed out. The MRR has a strong de- 
pendence on the number of points used to interpolate the 
values of the gravitational potential and its derivatives from 
the Fourier-grid on which they are calculated, as shown in 



the bottom row of Fig. (5] The use of more interpolation 
points allows for more robust ray-tracing, with a more ac- 
curate mean magnification value. However, the computation 
time rises approximately linearly with the number of inter- 
polation points used, and the interpolation scheme suffers 
in the vicinity of density peaks. The latter issue should not 
pose too much of a problem when the lensing mass distri- 
bution is large-scale structure, but computational efficiency 
presents a significant hurdle for ray-tracing procedures such 
as this. The parameters that most govern the accuracy of 
the ray-tracing scheme are the Fourier-grid size and mass 
smoothing length, which have essentially equivalent effects; 
assigning mass to an Fourier-grid acts to spread it out over 
a fixed number of grid-points. A CIC mass assigning scheme 
is different to a Gaussian filter smoothing, but in essence, 
the larger the area over which mass is spread, the more den- 
sity peaks are suppressed and the larger the error in the 
tests shown in Fig. [S] On the top row of Fig. [4l we show 
that halving the Runge-Kutta step-size has a negligible ef- 
fect, so half the Fourier-grid resolution is deemed sufficient. 
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Figure 2. The radial dependence of magnification {fi; left) and the percentage error in the numerical magnification (A/i//i; right) when 
ray-tracing through a discretized NFW lens. Different values of numerical parameters are compared. (Top) Radius of ray-bundle/image: 
9 = 0.1" (red circles) and 8 = 1" (blue squares); (Bottom) Number of points used in the interpolation scheme: ni„t = 2 (red circles) and 
^int = 1 (blue squares); In each case, 10* lines of sight have been distributed and binned uniformly with log(r/rvir), and Icr error-bars 
are shown. The analytic solution (black dots) at each impact parameter sampled is included in the panels on the left. In the panels on 
the right, the dotted horizontal line marks the 3 per cent error threshold. 



The model lens used for this test is subject to a few user- 
defined variables; a few convergence tests were run to ensure 
that these were not the cause of the observed features. Tri- 
alhng different mass resolutions from 10^ to W^^Mq for a 
fixed halo virial mass, M = lO^Mo, the results are shown 
in the panels of the middle row of Fig. 3] We find that the 
resulting scatter increases from less than 1 per cent up to 
10 per cent at a fixed impact parameter. One could equiva- 
lently say that the MRR is larger for haloes that are not as 
well resolved. The cosmological mass distributions in Sec. 15. II 
have mass resolutions better than lO^^Mo, so one can ex- 
pect a low error in magnification even within dense regions. 
The number of bins that the lens is divided into is a choice 
that has negligible effect on the ray-tracing results, as shown 
in the bottom row of Fig. 2] One is therefore reassured 
that the results of the ray-tracing on the discretised NFW 
profile are relevant to the results of ray-tracing through 
cosmological ma ss distributions based o n the fiducial pa- 
rameter choices. iKling fc Frittellil (|2008l ) numerically inte- 
grate the null geodesic equations (with the adaptive step- 
size Runge-Kutta-Fehlberg 4-5 method) in order to test the 



accuracy of the thin-lens approximation, but only with ref- 
erence to strong lensing by singul ar isothermal sphere (SIS) 
and Navarro- Frenk- White (NFW: lNavarro et al.lll995l 'l mass 
profiles, both of which are relatively thin. 



5 LARGE-SCALE STRUCTURE 
5.1 Cosmological simulations 

In order to construct the predicted probability distributions 
for the weak lensing statistics numerically, numerous null 
geodesic equations are integrated through a cosmological 
mass distribution that is generated with N-body simula- 
tions. The cosmological simulations are c arried ou t with 
the parallel Tree-PM SPH code GADGET2 l|Springell lioosh 
using collisionless particles only. We adopt a ACDM cos- 
mology, with the following values for the cosmological pa- 
rameters: f^M.o = 0.27, D.A,o = 0.73, /i=0.71, as = 0.9. The 
dark matter distribution is discretised into 256"^ particles 
distributed within a periodic box with co-moving length 
of I/box ~ 50 h~'^Mpc, resulting in a mass resolution of 
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Figure 3. The radial dependence of magnification (fi; left) and the percentage error in the numerical magnification [Afi/fi; right) when 
ray-tracing through a discretized NFW lens. Different values of numerical parameters are compared. (Top) Fourier-grid resolution: 24 
h~-'^kpc (red circles) and 49 h~^kpc (blue squares). (Bottom) No smoothing (red circles) and smoothing with a filter of width a = 1 
h~^kpc (blue squares). In each case, 10* lines of sight have been distributed and binned uniformly with log(r/rvir), and la error-bars 
are shown. The analytic solution (black dots) at each impact parameter sampled is included in the panels on the left. In the panels on 
the right, the dotted horizontal line marks the 3 per cent error threshold. 



TUp = 6.3 X 10 M©; the initial displacements are in a 'glass' 
configuration. The simulations begin at an initial redshift of 
Zi = 39, with a redshift-dependent gravitational force soft- 
ening length of Eco ~ 16 h~^kpc (Plummer-equivalent). 

The space between an observer and source is divided 
up into individual regions, each modelled by a snapshot of 
a cosmological simulation at an appropriate redshift. The 
snapshot cadence, Az, is chosen such that the hght travel 
time corresponds to the length of the boxes. The line-of- 
sight integrated co-moving distance between an observer and 
source at z = Zs is: 



dz 

'tm' 



(33) 



Since this cannot be easily evaluated in general, we instead 
assume that the differences between the redshift of adjacent 
boxes are small, and thus: 



cAz 

wry 



(34) 



which is easily rearranged to determine the snapshot ca- 



dence. A total of 28 boxes are required to fill the space be- 
tween an observer and a source at z « 0.5. 

Note that the scale factor at the position of a photon is 
required to integrate the null geodesic equations. This scale 
factor is not identified with the 'average' redshift of the cur- 
rent snapshot being traced through; instead it is derived 
from the lookback time (see Eqn. 13. 5[) and evolves during 
the integration through the box. That is to say, it is as- 
sumed that in the time it takes for a photon to traverse a 
box, the co-moving scale of the structure remains constant, 
however the physical scale length changes. As a sanity check, 
the numerically derived value of the scale factor at each box 
interface is compared to the redshift of the appropriate snap- 
shot. 



5.2 The sampling region 

In order to avoid repeated structure, a number of precau- 
tions were taken. Nine independent simulations were run, 
resulting in nine separate realisations. The snapshot for each 
required redshift was chosen randomly from among the nine 
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Figure 4. The radial dependence of magnification (p; left) and the percentage error in the numerical magnification [Afi/fi; right) when 
ray-tracing through a discretized NFW lens. Different values of numerical parameters are compared. (Top) Runge-Kutta step-size that 
is: 1/2 the Fourier-grid resolution (red circles) and for 1/4 the grid resolution (blue squares); (Middle) Mass resolution: rrip/MQ = 10^*^ 
(red circles), mp/M© = 10^^ (pink triangles) and rrip/M© = 10^ (blue squares); (Bottom) Number of radial bins used for the lens 
discretisation: 30 (red circles) and 10 (blue squares). In each case, 10'^ lines of sight have been distributed and binned uniformly with 
log(r/rvir), and la error-bars are shown. The analytic solution (black dots) at each impact parameter sampled is included in the panels 
on the left. In the panels on the right, the dotted horizontal line marks the 3 per cent error threshold. 
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Figure 5. A two-dimensional schematic diagram of the ray-tracing method. Large simulation snapshots (dashed lines; three shown here) 
of length -Lbox fiU up the space from the observer, 2 = 0, to the source plane at Zs- Small cubes (thick lines) from within these boxes 
form a prism of lensing mass. The mass in each cube is distributed onto Fourier-grids (thin lines) before calculating the gravitational 
potential and derivatives in order to numerically integrate the path (green) of a photon. 



realisations. The chosen box was then randomly translated 
and rotated 90° about any /all of the axes. From each 50 
h~^Mpc box we selected the lensing mass formed by par- 
ticles within a prism 6.25 h~^Mpc across the sky and 50 
h~^Mpc along the line of sight. This is split into 8 cubes 
each of volume V = 6.25^ (h-iMpc)^ ^ 244 (h~^Mpc)^ 
Each of these has its mass placed on a 128'^ point grid, us- 
ing the CIC algorithm. We apply an FFT to these boxes 
with zero-padding so that the density field (and the deriva- 
tives of the gravitational potential) due to structure within 
the box only is determined on a scale of 49 h~^kpc. Note 
here that the simulation is run on a larger box such that the 
large-scale modes are included in the formation of structure, 
but we only use a smaller portion of the box for the lensing 
(see Fig. [5] for a cartoon diagram). The null geodesic equa- 
tions of the eight rays and one anchor of each ray-bundle 
are numerically integrated using the methods described in 
Sec. 13.41 Lines of sight that fall too close to the edge of the 
boxes will be artificially sheared, so assuming that most of 
the lensing occurs due to mass within 0.5 h^^Mpc, we only 
send out 'beams' that will fall within the central 5.25 x 5.25 
(h~^Mpc)^ region at the source redshift. For a source red- 
shift of z=0.5, with an angular diameter distance of approxi- 
mately 900 h~^Mpc, this means that the area of sky sampled 
is about 0.33° x 0.33° or 5 x 10"" of the entire sky. To in- 
crease the sampling area, this entire procedure is used to 
trace the paths of 10 000 bundles, then repeated for another 
10 000 bundles using another lensing mass distribution sam- 
pled from the same set of simulations. 

5.3 Numerically determining the /^PDFs 

Using ray-tracing procedures to construct a magnification 
probability distribution is effectively equivalent to calculat- 
ing the angular diameter distance in different directions out 
to the same co-moving distance. The distribution of angu- 



lar diameter di stances for a given source redshift has been 
investigated bv lHadrovic fc BinnevI (|l997h . In the extremal 
case of a large beam-size, and therefore, large separations 
between rays, the angular diameter distances would be the 
same in all directions and the 'distribution' would be a 
single peak at Dfb{z), the average FLRW value, implying 
homogeneity on those scales. However, if the beam-size is 
small enough, then the inhomogeneities induce a distribu- 
tion. Each ray-bundle probes a single light of sight. More 
specifically, it describes the lensing that has affected a sin- 
gle image of a fixed size on the sky. The number of bun- 
dles that exhibit a magnification of fi are thus proportional 
to the probability that an image is magnified by jj,. How- 
ever, the statistical quantity of interest is the magnification 
probability distribution for sources. Therefore, the number 
counts that are used to produce the magnification probabil- 
ity histograms presented in Sec. l5.4l are weighted by the area 
covered by the beam at the source plane, or equivalently, the 
inverse of the magnification. 

where F(/ii,^2) is the fraction of all ray-bundles for which 
Hi < fj. < jj,2- Similarly the mean and standard deviation of 
these properties are also weighted. 

5.4 Comparison of the PDFs produced with and 
without lens planes 

The predicted probability distributions of magnification and 
shear are constructed by ray-tracing 2 x 10* bundles through 
the three-dimensional mass distribution described above. 
The results, shown in Figures |6] and [T] are compared to 
the probability distributions constructed with the multiple 
lens-plane RBM using 5 x 10* bundles. However, not all ray- 
bundles are included in the analysis. The RBM is very well 
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Figure 6. The magnification probabiUty distribution for a source 
at redshift z = 0.5 as predicted by the multiple lens-plane RBM 
(blue dashed line) and by the three-dimensional equivalent (black 
solid line) 



Figure 7. The shear probability distribution for a source at red- 
shift z = 0.5 as predicted by the multiple lens-plane RBM (blue 
dashed line) and by the three-dimensional equivalent (black solid 
line) 



suited to the weak lensing limit and provides a computation- 
ally efficient alternative to grid-based ray-shooting methods; 
however, it underestimates the magnification in cases where 
multiple imaging is expected as only one image contributes 
to the magnification of each source. For this reason, rays that 
fall within a minimum distance of a lens (remembering that 
a distribution of point masses describes the lensing mass) 
are excluded from the analysis. This affects the high magni- 
fication probabilities deduced, but does not sign ificantly af- 
fect the weak lensing analysis (jFluke et al.ll2002l ). The RBM 
avoids the danger of artificial shear due to a source pixel col- 
lecting light rays that have passed near the edge of a shoot- 
ing grid. If a ray-bundle were to suffer from neither conver- 
gence nor shear, it would obtain the minimal magnification, 
/ie(,,min = 1, Or equivaleutly, fift^^in = D%{zs)/Dl,,{zs) (see 
AppendixlEj. For a source redshift of 2 = 0.5, this minimum 
magnification is ~ 0.965; the relevant angular diameter dis- 
tances h ave been determin ed with the aid of the Angsiz 
routine (|Kavser et all ligPTh . Those bundles that produce 
magnifications below this minimum value represent the de- 
magnified components of a multiple-image system, a possible 
consequence of strong lensing. These will highly underesti- 
mate the total magnification of the associated source, and 
so, are also excluded in both methods. On the other hand, 
if a magnified image is traced back to the source plane, but 
belongs to a multiple-image system, it will not be identi- 
fied as such. Although it will also underestimate the total 
magnification, it will generally do so by a negligible amount. 
The ray-tracing method developed in the present work does 
not produce any bundles that must be excluded. The av- 
erage magnification predicted by both methods is close to 
unity, as required by fiux conservation (see Egn. I14p . More 
precisely, the multiple-lens plane approach finds fl = 1.0199 
while the three-dimensional approach finds a mean magnifi- 
cation of p. = 0.9994. Using the three-dimensional method to 
analyse lensing by large-scale structure, we find that large 
magnifications are not produced, as shown in Fig. [B] The 
slope of the differential magnification probability distribu- 
tion is found to be much steeper than for RBM for the regime 



1.1 < ^ < 1.5. The difference in the two methods can also 
quantified by the standard deviation, found to be — 0.205 
from the multiple-lens plane approach and ct^ — 0.021 from 
the three-dimensional approach, and order of magnitude 
lower. To understand the reason of this difference, we also 
plot the probability distribution for shear, shown in Fig. [71 
While the overall shape of the distribution function is sim- 
ilar as derived by the two methods, the multiple-lens plane 
method measures more high shear values and less low shear 
values. The projection of matter onto multiple lens-planes 
may have artificially increased the shearing effect of cosmo- 
logical lenses. However, the effect of shot-noise in the stan- 
dard RBM method is identified as another possible reason 
for the discrepancy. When the three-dimensional approach 
is used, the gravitational potential is calculated by means 
of Fourier-methods. In contrast, the multiple-lens plane ap- 
proach determines the potential field at the location of each 
ray, but adding the potential due to each simulation parti- 
cle in the plane individually. Each particle in the simulation 
represents some underlying smooth mass distribution, but 
by treating them as point masses, t he method is susc e ptible 
to the effects of shot-noise. Recentlv. lTakahashi et al] (|201ll ) 
discussed the effects of shot-noise on the variance of the con- 
vergence found by ray-tracing through N-body simulations. 
Although they did not use the same ray-bundle approach 
that we do here, they apply a multiple lens-plane method 
with an FFT for calculating the gravitational potential in 
each plane; they compare their results for a range of Fourier- 
grid resolutions (see Figure 2 of their paper). When the grid 
resolution is large enough to smear out structure, the vari- 
ance of the convergence falls below theoretical predictions; 
although our choice of statistic is different, our results agree 
with this interpretation (see Fig.[3]in Sec. 14. 21) . Interestingly, 
a very small grid-resolution (< 5 h^^kpc) leads to a variance 
that is larger than the theoretical prediction, a result which 
they attribute to shot-noise. We agree with this interpreta- 
tion, noting that Fourier methods applied to grid-resolutions 
that are much smaller than the mean-interparticle distance 
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Figure 8. The magnification probabiUty distribution for a source 
at redshift z = 0.5 as predicted by the three-dimensional method. 
Two different force-softening lengths applied in the N-body simu- 
lations are compared: e = 5 h~^kpc (red dashed line), and e = 16 
h~^kpc (black solid line) 



would resemble a direct summation as employed by the 
multiple-lens plane approach in the present work. 

The /iPDF predicted by the three-dimensional method 
is subject to certain choices of numerical parameters, some 
of which have been discussed above. Thus, we test the pa- 
rameters that are most likely to have an impact on the cos- 
mological lensing results presented above. For example, a 
parameter that was not relevant to the tests presented in 
Sec. m is the force-softening length chosen for the cosmolog- 
ical N-body simulation. We re-run the simulation described 
earlier but with a force-softening length reduced to e = 5 
h~^kpc. We sample 3 x 10* lines of sight through the mass 
distribution determined with this simulation, and in Fig. |8l 
we compare the /iPDF to the result for the three-dimensional 
method shown in Fig. [S] For each, the Fourier-grid resolu- 
tion is 49 h~^kpc. The similarities in the PDFs reassures 
us that any smoothing of structure below 16 h~^kpc is not 
responsible for any error in the PDFs measured by the three- 
dimensional method. 

In addition, as the Fourier-grid scale is reduced, the 
numerical solution is expected to be more accurate, with the 
caveat that it remains large enough to prevent the effects of 
shot-noise. In Fig. [5] we compare two different Fourier-grid 
resolutions: 24 h~^kpc for which 9 x 10* lines of sight have 
been sampled, and 49 h~^kpc for which 3 x 10* lines of sight 
have been sampled. Both tests are run on simulations with 
a force-softening length reduced to e = 5 h~^kpc, which are 
only able to reliably describe structure on scales larger than 
20 h~^kpc. The differences are negligible, which demonstrate 
that most of the weak lensing results from structure above 
49 h~^kpc scales. Structures on scales larger than 24 h~'^kpc 
that are able to produce strong or even 'moderate' lensing 
are rare. 

We recognize that there are multiple ways to create the 
same Fourier-grid resolution. If the number of cubes taken 
from the simulation boxes are doubled, but the number of 
Fourier grid-points halved, the resolution does not change. 
However, only a quarter of the patch of sky would be sam- 
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Figure 9. The magnification probability distribution for a source 
at redshift z = 0.5 as predicted by the three-dimensional method. 
Two different Fourier-grid resolutions are compared: 24 h~^kpc 
(green solid line), and 49 h~'^kpc (red dashed line) 
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Figure 10. The magnification probability distribution for a 
source at redshift z = 0.5 as predicted by the three-dimensional 
method. The same Fourier-grid resolution (24 h~^kpc) is created, 
but excising either 8 cubes (purple dashed line) or 16 cubes (yel- 
low solid line) from each simulation box 



pled, and the computational run-time would double. On the 
upside, the memory usage reduces to an eighth, which is a 
significant advantage for such a computational demanding 
approach. In Fig.[TO]we show the result of doubling the num- 
ber of cubes in such a manner; the Fourier-grid resolution 
is 24 h~ kpc. A total of 3 x 10 lines of sight have been 
sampled for the method using 8 cubes, while 6 x 10* lines 
of sight have been sampled for the method using 16 cubes. 
Note that this was not tested on the compact lens models in 
Sec. |4] since cubes are not excised from simulations as they 
are when analysing cosmological lensing. We are satisfied 
that the choice of 8 cubes for our previous results (see Fig. [6] 
and Fig. [7)l has no bearing on the results. 
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6 SUMMARY AND DISCUSSION 

Modelling the magnification probability distribution of 
background sources due to gravitational lensing relies on 
the application of ray-tracing methods. This work has taken 
the first steps towards making a direct comparison between 
the predictions of the multiple lens-plane ray-bundle method 
(RBM) and one that does not invoke the thin-lens approxi- 
mation; instead, the null geodesic equations are integrated. 
The efficiency and accuracy of this computationally chal- 
lenging approach can be improved by careful choices of nu- 
merical parameters; therefore, the results are analysed for 
the behaviour of the ray-tracing code in the vicinity of 
Schwarzschild and NFW lenses. A range of tests were able to 
pin down the numerical parameters that play a critical role 
in the predicted statistics. The behaviour of the ray-bundles 
in the vicinity of a Schwarzschild lens demonstrated that the 
method can reproduce large magnifications given sufficient 
spatial resolution. The limitations are dominated by the spa- 
tial resolution of the Fourier-grid and the mass resolution of 
the discretized lens. 

Comparisons to a multiple lens-plane algorithm are 
drawn in the context of cosmological mass distribution for 
a source redshift of Zs = 0.5. The weak lensing statis- 
tics predicted by ray-tracing through simulated cosmologi- 
cal mass distributions found significant differences compared 
to the results from the original RBM approach. Either the 
use of multiple lens-planes or shot-noise are responsible for 
the observed differences. To clarify the dominant factor, a 
multiple lens-plane ray-bundle method that applies a two- 
dimensional Fourier-method at each lens-plane is proposed 
as the next step towards the direct comparison. 

The method developed here presents a computational 
challenge as the three-dimensional FFTs prove memory- 
expensive, but is justified by the need to quantify the mul- 
tiple lens-plane approximation. We are, as yet, unable to 
sample the large number of lines of sight required to model 
the intermed iate magnification {2 < fi < 20) region where 
iRauch et al.l ( jf 992 ) had identified the feature present only 
in two dimensional lens models. 

We now discuss future applications of the this ray- 
tracing method that are outside the scope of the present 
work. For example, proposed large surveys have the mea- 
surement of cosmic shear, with particular focus on the de- 
termination of the nature of dark energy, as one of their main 
science drivers, using a combination of a large area of sky 
coverage and high-precision photometric redshifts. These in- 
clude the ground-based Panoramic Survey Telescope and 
Rapid Response System (Pan-STARRSfl the VST-KIlo- 
Degree Survey (KIDSfl the Dark Energy Survey (DESfl, 
the Large Synoptic Sjrrvoy Telescope (LSST jf| as well as the 
space- based Eucfid iLaurciis 2009f| and the Joint Dark En- 
ergy Mission (JDEM)Q, which has recently been rebranded 
as the Wide- Field InfraRed Survey Telescope (WFIRST). 
The method developed in the present work has the potential 



http:/ /pan-Starrs. ifa.hawaii.edu 
^ http:/ /www. astro-wise.org/projects/KIDS/ 
■* http:/ /www. darkenergysurvey.org/ 
^ http://www.lsst.org 
^ http://sci.esa.int/euclid 

^ http://jdem.lbl.gov/ and |http: //jdem.gsfc.nasa.gov/| 



to quantify the effects of multiple lens-plane techniques on 
the theoretical predictions for cosmic shear measurements. 
Another alternative to ray-tracing is to retain the three di- 
mensional mass distribution, but assume that the defiections 
are small enough to satisfy the Born approximation. Here, 
the matter that is integrated along the (un-defiected) line of 
sight is solely responsible for the convergence. Over larger 
and larger distances, the Born approximation becomes less 
and less accurate; several studies have found the requirement 
for corrections in the construction of t he shear power spec- 
trum and higher order bi spcctrum (e.g. lVan Waerbeke et all 
I2OOII : IShapiro fc Cooravl 2006). We envisage that the three- 
dimensional method could be used to make a one-to-one 
comparison with the Born approximation, just as we have 
done here with the multiple lens-plane method. 

Finally, with a larger number of ray-bundles and greater 
sky-sampling, we could turn our attention to higher-order 
statistics. Flexion is the third-order effect in gravitational 
lens ing and is effectively the g radient of the shear component 
(see IColdberg fc BaconlboOSh . It has recently been demon- 
strated that flexion ca n be modelled using a ray-bundle 
method (Fluke fc Laskvll201l| j: this presents an exciting fu- 
ture application of the method, but surely it is even more 
necessary to justify the use of the multiple lens-plane for- 
malism. 
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APPENDIX A: CHRISTOFFEL SYMBOLS 

The expanding weak-field metric presented in Eqn. [18] has 
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Note that x denotes , y denotes x^, and z denotes x'^. 
The ChristofFel symbols for this metric are dependent, not 
only on the gradients of the gravitational potential <p, but 
(j) itself. However, the weak field condition requires that the 
magnitude of the perturbations satisfy <C a. The Geodesic 
Equations - the second order differential equations for the 
four coordinates - are then constructed. 



the lens, so the convergence is zero. The solution for the 
shear for images lensed by a Schwarzschild lens, 7schw, is: 

7Schw(2:) = X. (C3) 

The analytic solution for the magnification for an image at 
any location, /ischw, is: 




APPENDIX B: 
SCHEME 



THE INTERPOLATION 



Here we only describe the scheme for the evaluation of 
(f>{xl,x^,x^) given the gridded values (f){x^ , x^ , x^) , but the 
same method is applied to evaluate the derivatives of the 
potential. The scheme is not strictly tri-cubic interpolation, 
but a three-dimensional version of the bicubic spline. In- 
stead, it initially performs a one-dimensional spline to com- 
pute the derivative with respect to the line of sight; i.e. it 
fits a spline to find d(l){x^ , x^ , x"^) /dx'^ across the entire grid. 
This is only done once for each lensing mass distribution. 
Each time the local values need to be found, the scheme 
identifies the gridded values that include runt points on ei- 
ther side of the current location for each dimension, i.e. a 



steps: 



mesh grid. Then, the scheme performs the following 



(i) The gridded derivative is interpolated across the 
chosen dimension, in this case the a;^-dimension, to find 
4>{x^ , x^, xl), i.e. the value with the current x'^-coordinate of 
the photon specified, but the other coordinates correspond- 
ing to grid points. (2ni„t)^ interpolations are required. 

(ii) Another set of 2nint one-dimensional splines 
are performed, this time across the a;^-axis, to find 
d(t){x'^,x'^,xl)/dx^. 

(iii) 2ni„t interpolations along the x^-axis determine 
(t){xl,x^,x'i). 

(iv) A single one-dimensional spline across the x^-axis is 
used to find d4>{xl, x'^ , xl) /dx'^ 

(v) Finally, a single interpolation along the x^-axis allows 



APPENDIX C: SCHWARZSCHILD LENS 

A single point of infinite density is one of the simplest lens 
models, and is often used to model a spherically symmetric 
lens for idealised analytical studies. The relevant (angular) 
scale length for a so-called Schwarzschild lens of mass AI is 
the Einstein radius: 



D.Ds c 



(CI) 



At the lens plane, this corresponds to a physical scale length 
given by: 



(C2) 



The gravitational lensing quantities of interest are found by 
combining the derivatives of this potential, as described in 
Sec. (2] They are defined in terms of a dimensionless radial 
distance, x = r/vEin- All lines of sight probe regions outside 



APPENDIX D: NFW LENS 

The NFW profile, given in Eqn. 1311 can be considered a 
thin lens for cases where the observer, lens and source are 
separated by large angular diameter distances. The follow- 
ing analytic solutions, ID1IID81 are therefore derived from 
projected surface densities. A dimensionless radial distance, 
X = r/ra, has been adopted. Here, the quantities pc, r^, Sc 
an d Serif are those defined in Eqns. l22|[3l]l32l andl4l Follow- 
ing l^ighTfeBrainer^ |200d ) and lCoi (iSjloi ). the analytic 
solution for the radial dependence of the convergence is given 
by 



I^Y^¥W{X) = — K{x) 



(Dl) 



where K{x) is given by 
1 - 



K{x) = <^ 
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■■ tan 



if a; < 1 
if a; = 1 
iix>l 



By integrating Eqn. lDll over the area within r, one finds the 
mass within a cylinder of radius r: 



MNPW,cyl(a;) = 'iTTT-aScPcCix), 

where C{x) is given by 
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The radial dependence of shear is thus given by 
7NFw(a;) = (j\x), 



where Gix) is given by 

9<'\x) 



G(x) 



if a; < 1 
if a: = 1 
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gy (x) if a; > 1 

The functions g<c{x) and g>{x) are given by 
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and 
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Combining these with Egns. iDll and lDSl with Eqns. [3]and[8] 
allows one to determine the analytic value of the magnifica- 
tion, /X, of an image at a given projected distance from the 
centre of the NFW lens. 



APPENDIX E: FULL BEAM VS EMPTY BEAM 

Magnification is the relative increase in flux as a result of 
gravitational lensing. There is a subtlety here. One may com- 
pare the flux received from the source to the flux received 
if the universe was entirely homogeneous, in which case the 
magnification is referred to as the full-beam magnification, 
/x/i,. This is the magnification referred to in Sec. 12.11 The 
so-called empty-beam magnification, jj,eb, is defined relative 
to the case where all the matter in the universe is locked up 
in compact objects, and all lines of sight are empty. As this 
is the case of minimal magnification (no convergence and 
assumed negligible shear), then n^t is always greater than 
unity. To derive the relationship between the two magnifi- 
cations, we must compare the solid angles subtended by the 
source: ilfb in the full-beam scenario, il^b in the empty-beam 
scenario and Qimg the solid angle of the image observed. 
Given the conservation of surface brightness, the magnifica- 
tions are defined by: 

fifb = — — and = -jT-^- (El) 

Si/6 i'eb 

The solid angles are related to the physical area of the source 
via the angular diameter distances: 

(1 + zyDfb(z) 

is the appropriate relationship for the full-beam case, and 

"-(-^- (l + .m.(.) - 

is appropriate for the empty-beam case. Combining Equa- 
tions lEll IE2I and IE31 one is easily able to convert between 
the two: 

D^^,{z)- ^^^^ 

The angular diameter distance in the full-beam case, Dfb{z), 
is equivalent to the Dyer-Roeder distance D{a = 1; z) which 
is also the FLRW solution: 

c dz' 
^''^'^^ Ho{l + z)j, [fi^,„(i + ,,)3 + j^^^]i/2- (E5) 

The empty-beam angular diameter distance, Deb{z), is 
found by solving the Dyer-Roeder equation for 3 = 0: 



Deb{z)^D(&^Q-z) = cj^ (TT^T^T^rf-' 



(E6) 



This paper has been typeset from a T^jX/ I^T^^X file prepared 
by the author. 



